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Abstract 

We study a Fokker-Planck equation modelling the firing rates of 
two interacting populations of neurons. This model arises in com- 
putational neuroscience when considering, for example, bistable visual 
perception problems and is based on a stochastic Wilson-Cowan system 
of differential equations. In a previous work, the slow-fast behavior of 
the solution of the Fokker-Planck equation has been highlighted. Our 
aim is to demonstrate that the complexity of the model can be dras- 
tically reduced using this slow-fast structure. In fact, we can derive a 
one-dimensional Fokker-Planck equation that describes the evolution 
of the solution along the so-called slow manifold. This permits to have 
a direct efficient determination of the equilibrium state and its effec- 
tive potential, and thus to investigate its dependencies with respect to 
various parameters of the model. It also allows to obtain information 
about the time escaping behavior. The results obtained for the reduced 
ID equation are validated with those of the original 2D equation both 
for equilibrium and transient behavior. 

1 Introduction 

In this work, we will propose a procedure to reduce rate models for neuron 
dynamics to effective one dimensional Fokker-Planck equations. These sim- 
plified descriptions will be obtained using the structure of the underlying 
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stochastic dynamical system. We will emphasize the numerical and practi- 
cal performance of this procedure coming from ideas used in the probability 
community [''>] for a particular model widely studied in the computational 
neuroscience literature. 

We will consider a simple model [15, 13, 12] formed by two interacting 
families of neurons. We assume that there is a recurrent excitation with a 
higher correlation to the activity of those neurons of the same family than 
those of the other while a global inhibition on the whole ensemble is due 
to the background activity. These families of neurons are modelled through 
the dynamics of their firing rate equations as in the classical Wilson-Cowan 
equations [J 7]. The synaptic connection coefficients Wij, representing the 
strength of the interaction between population i and j, are the elements of 
a 2 X 2 symmetric matrix W given by 



W 



W+ — Wj W- — Wj 

W- — Wl w+ — Wj 



Here, Wj^ is the self-excitation of each family, W- the excitation produced on 
the other family, and wi the strength of the global inhibition. The typical 
synaptic values considered in these works are such that W- < wj < 
leading to cross-inhibition since W- < wj and self-excitation since wj < 
Let us comment that these rate models are very simplified descriptions 
of interacting neuron pools, more accurate microscopic models introducing 
neuron descriptions at the level of voltage and/or conductances probability 
distribution can be derived [6, 7, 8, 9, 10, 11]. 

The time evolution for the firing rates Ui{t) of the neuronal populations 
i = 1, 2 as given in [' 1] follows the following Stochastic Differential Equations 
(SDE): 



= -v^{t) + (t>[Xi+Y, WijUjit) ] + PMt), i = 1, 2, (1) 



dt 



where r = 10~^s is a time relaxation coefficient, which will be chosen equal 
to 1 in the sequel except for the numerical results, and (,i{t) represents a 
white noise of normalized standard deviation Pi > 0. In (1) the function 
(f){x) is a sigmoid function determining the response function of the neuron 
population to a mean excitation x{t) = Xi + Ylj ^ij'^j- 

cPix) 



1 -|- exp(— a(x/fc — 1)) 
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where Aj are the external stimuh apphed to each neuron population. 

We will recall in the next section that the study of the decision making 
process for the previous network can be alternatively studied by means of 
the evolution of a Fokker-Planck equation in two dimensions i.e. the plane 
{t^i-,i^2)- The theoretical study of such problem (existence and uniqueness 
of positive solutions) was done in [o]. However, we will emphasize that due 
to slow-fast character of the underlying dynamical system the convergence 
towards the stationary state for the corresponding two-dimensional problem 
is very slow leading to a kind of metastable behavior for the transients. 

Nevertheless, the 2D Fokker-Planck equation allows us to compute real 
transients of the network showing this metastable behavior. Moreover, we 
can derive a simplified one dimensional SDE in Section 3 by scaling con- 
veniently the variables. Here, we use the spectral decomposition and the 
linearized slow manifold associated to some stable/unstable fixed point of 
the deterministic dynamical system. The obtained ID Fokker-Planck equa- 
tion leads to a simple problem to solve both theoretically for the stationary 
states and numerically for the transients. In this manner, we can reduce the 
dynamics on the slow manifold to the movement of a particle in an effec- 
tive ID potential with noise. We recover the slow-fast behavior in this ID 
reduction but, due to dimension, we can efficiently compute its numerical 
solution for much larger times than the 2D. We can also directly compute an 
approximation to the 2D equilibrium state by the ID equilibrium onto the 
slow manifold since in ID, every drift derives from a potential. 

Let us mention that another approach to get an approximation of the 2D 
Fokker-Planck equation by ID Fokker-Planck reduced dynamics has been 
proposed in [16]. This approach is purely local via Taylor expansion around 
the bifurcation point leading to a cubic ID effective potential. Moreover, an 
assumption about the scaling of the noise term is performed to be able to 
close the expansion around the bifurcation point. Our approach is valid no 
matter how far we are from the bifurcation point as long as the system has 
the slow-fast character and we do not assume any knowledge of the scaling 
of the noise term. Moreover, we can reconstruct the full potential not only 
locally at the bifurcation point. We point out that the results of their ID 
Fokker-Planck reduction are compared to experimental data in [Hi] with 
excellent results near the bifurcation point. A similar applied analysis of our 
reduced Fokker-Planck dynamics in a system of interest in computational 
nuroscience is underway [ .]. 

Section 4 is devoted to show comparisons between the 2D and the reduced 
ID Fokker-Planck equations both for the stationary states and the transients. 
We demonstrate the power of this ID reduction in the comparison between 
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projected marginals on each firing rate variable and on the slow linearized 
manifold. Finally, section 5 is devoted to obtain information of the simulation 
in terms of escaping times from a decision state and performance in the 
decision taken. 

2 The two dimensional model 

We will illustrate all our results by numerical simulations performed with the 
physiological values introduced in [j2]: a = 4 and Vc = 20Hz, Ai = 15Hz 
and A2 = Ai+AA, with AA = for the unbiased case and AA = 0.01, 0.05, 0.1 
for the biased case. The noise parameter is chosen as f3 = 0.1, and the 
connection coefficients are given by w+ = 2.35, wj = 1.9 and W- = 1 — 
r{w-^ — 1)/(1 — r) with r = 0.3. 

It is well known [12] that the deterministic dynamical system associated 
with (1) is characterized by a supercritical pitchfork bifurcation in terms 
of the parameter w+ from a single stable asymptotic state to a two stable 
and one unstable equilibrium points. We recall that the unstable point is 
usually called spontaneous state while the two asymptotically stable points 
are called decision states. The behavior of the bifurcation diagram for the 
deterministic dynamical system defining the equilibrium points in terms of 
the Wj^ parameter and with respect to AA is shown in Figure 1. Observe 
that in the nonsymmetric (AA 7^ 0) bifurcations, the pair of stable/unstable 
equilibrium points detaches from the branch of stable points. 



Figure 1: Bifurcation diagram: z/i-component of the equilibrium states with 
respect to w+. Left Figure: bifurcation diagram for the unbiased case AA = 
0. Right Figure: bifurcation diagram for the biased case AA = 0.1. 

For example, with w+ = 2.35 in the unbiased case, if AA = 0, the stable 
points are in 5*1 = (1.32,5.97) and its symmetric S3 = (5.97,1.32), and the 
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unstable one is in S2 = (3.19, 3.19); whereas, in the biased case AA = 0.1 the 
stable points are in Si = (1.09,6.59) and ^3 = (5.57, 1.53), and the unstable 
one in S2 = (3.49,3.08). 

Furthermore, it can be shown by means of direct simulations of system 
(1), that there is a slow-fast behaviour of the solutions toward the equilibrium 
points. This behavior is plotted in Figure 2, where the straight lines show 
the behavior of several realisations for the deterministic system (i.e. when 
Pi = 0), and the wiggled line represent one realisation for the full stochastic 
system (1). Figure 2 highlights also the so called slow manifold: a curves in 
which the three equilibrium points of the system lie and where the dynamics 
are reduced to rather quickly. 




Figure 2: Dynamics of a firing rate towards stable equilibrium, fast conver- 
gence towards the slow manifold and slow convergence towards one of the 
stable equilibrium points along the slow manifold. 

Applying standard methods of Ito calculus, see [14], we can prove that 
the probability density p = p(t, z^), with t > and = {i'i,i'2), satisfies 
a Fokker-Planck equation (known also as the progressive Kolmogorov equa- 
tion). Hence, p(i, i^) must satisfy: 

dtp + V-{{-u + <^>{A + W-i^))p)-^Ap = (2) 

where £ n = [0,1/^,] x [0, f^], A = (Ai,A2), $(xi,X2) = (0(xi), 0(x2)), 
V = {dui,du2) and A = Ay. We complete equation (2) by the following 
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Robin boundary conditions or no flux conditions: 

(-z/ + $)p- yVpj -71 = (3) 

where n is the outward unit normal to the domain. 

Physically, this kind of boundary conditions means that we have no par- 
ticles incoming in the domain. This is naturally relevant for the boundaries 
[0, I'm] X {0} and {0} x [0, Um]- For the two others boundaries, [0, Vm] x {vm} 
and {um} X [0, fm], it relies on the choice of large enough in such a way 
that the evolution of our system of particles is isolated. In practice, for the 
choosen parameters of our model, z/^ = 10 is a good choice. 

In order to simplify notations, let us consider, from now on, the vector 
field F = (Fi, F2) representing the flux in the Fokker-Planck equation: 



then, equation (2) and boundary conditions (3) read: 

9tP + V- (^Fp-^Vp^ =0 (5) 

Fp - yVp] - 71 = (6) 

We refer to [•]] for numerical results and a detailed mathematical anal- 
ysis of the Fokker-Planck model (5)-(6): proof of the existence, uniqueness, 
and positivity of the solution, and its exponential convergence towards the 
equilibrium, or stationary state. Let us just recall that the equilibrium state 
cannot be analytically given because the flux does not derive from a poten- 
tial, i.e. it is not in gradient form. 

Moreover, we remark that the slow-fast structure leads to stiff terms and 
thus, to small time steps and large computational time. In fact, the slow 
exponential decay to equilibrium makes impossible to wait for time evolving 
computations to reach the real equilibrium. Hence, it is difficult to nu- 
merically analyze the effect of the various parameters of the model on the 
equilibrium state, and then the importance of deriving a simplified model 
capable of explaining the main dynamics of the original one is justified. Nev- 
ertheless, one could find the equilibrium state directly by numerical methods 
to find eigenfunctions of elliptic equations. The discussed slow-fast behavior 
will serve us, in the sequel, to reduce the dynamics of the system to a one 
dimensional Fokker-Planck equation. 



l/l -I- (^(Ai -I- W\\V\ -I- W\2V2) \ 
V2 + ^{\2 + W2\V\ + W22l^2) J 
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3 One dimensional reduction 



In this section we present the one dimensional reduction of system (1). We 
shall treat first the deterministic part, see 3.1, then the stochastic terms, 
section 3.2, and finally we describe the one dimensional Fokker-Planck model, 
see 3.3. 

3.1 Deterministic dynamical system 

The slow-fast behavior can be characterized by considering the deterministic 
system of two ordinary differential equations, i.e. (1) with /3j = 0. Regardless 
of the stability character of the fixed point 5*2 = {vl'^ , 1^2'^), the slow-fast 
behavior is characterized by a large condition number for the Jacobian of 
the linearized system at the equilibrium point S2, i.e., a small ratio between 
the smallest and largest eigenvalue in amplitude. 

More precisely, let us write the deterministic part of the dynamical system 
(1) as follows: 



where z/ is a vector and F(l') = —v -\- $(A -|- Wv) is the flux, see (4), as 
described in section 2. Let us denote v^'^ the spontaneous equilibrium point, 
so that F{u'^'^) = 0. By spontaneous state we mean the only equilibrium 
before the bifurcation point and the unique unstable equilibrium point af- 
ter the subcritical pitchfork bifurcation. This equilibrium point v'^'^ is then 
parameterized by the bifurcation parameter Wj^ and it has a jump discon- 
tinuity at the bifurcation point for nonsymmetric cases A A 7^ 0. Hence, by 
construction: 



where we have denoted by Zi the values Zi = \i + wnvi + Wi2V2- 

We recall that v'^'^ is an hyperbolic fixed point (saddle point) after the 
bifurcation while before it is an asymptotically stable equilibrium. Hence 
the Jacobian Jp{v'^'^) has two real eigenvalues ^1 and ji2 being both neg- 
ative before the bifurcation and of opposite signs after. The bifurcation is 
characterized by the point in which the smallest in magnitude eigenvalue 
becomes zero. Let us denote by jii the (large) negative eigenvalue and by 
//2 the (small) negative/positive eigenvalue of Jf{v'^'^)- We remark that, the 



(7) 



{ul\ul^) = ^{K + W{ul\ul^)). 
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small parameter e << 1 which is responsible for the slow-fast behavior is 
determined by the ratio of the amplitude of the two eigenvalues: 



The values in terms of w+ and for different values of A A are shown in figure 
3. In the range of parameters we are interested with, e is of the order of 
10"2. Note the jump discontinuities at the bifurcation point for AA 7^ 
since the point around which our analysis can be performed jumps to the 
new created branch of the bifurcation diagram at the bifurcation point. 




Figure 3: Eigenvalues ratio with respect to W-^. and for AA = 0, 0.01, 0.05, 0.1. 

In order to reduce the system we need to introduce a new phase space 
based on the linearization of the problem. We will denote by P the matrix 
containing the normalized eigenvectors and by its inverse matrix. Note 
that, in the unbiased case (AA = 0), we have: 

with the associated eigenvalues fii = —1.55 and /X2 = 0.036, and the eigen- 
vectors are orthogonal. Orthogonality of the eigenvectors is no longer true 
for the nonsymmetric biased problem AA 7^ 0. Furthermore, using Hartman- 
Grobman theorem [J , 2], we know that the solutions of the dynamical system 
are topologically conjugate with its linearization in the neighbourhood of an 
hyperbolic fixed point, which is valid in our case for all values of the bi- 
furcation parameter except at the pitchfork bifurcation. Let us write it as 
follows: 

Jf{u^1) = PDP-\ (10) 
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where P is the matrix of eigenvectors and D is the associated diagonal ma- 
trix. We can describe the coordinates v in the eigenvector basis and centered 
on the saddle point u^'^ as follows: 

v = v'"i + PX, (11) 

which gives the definition for the new variable X = see also figure 4: 




Figure 4: Change of variable from the phase space {vi,V2) to {x,y). 

We can conclude that system (7) reads in the X phase space as: 

X = H{X) (12) 

where H{X) is the two dimensional vector defined by : 

H{X) = p-^Fiu^" + PX). 

We remark that by means of the chain rule, the Jacobian Jh{X) is given by: 

Jh{X) = p-^JFii^"" + PX)P, 

and using (10) and that X{i''^'^) = 0, we obtain that Ji^(O) = D, which is 
the diagonal matrix in the change of variables (10). 

Let us now make explicit the system (12) in terms of its components 
f{x,y) and g{x,y): 

X = f{x,y) 

y = 9{x,y) 
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where considering the definition of the flux F given by (4): 



fix, y) = -x- {p-^v'^)i + (P-1«>(A + W{v''i + PX)))i (13) 

g[x, y) = -y- (P" V^^)2 + [p-^^h + Wiv'"^ + PX)))^ 

Now, we can choose a new time scale for the fast variable r = et, with e 
given by (8), in such a way that for large time t ~ then the fast variable 
T ~ 1 and the variations ^ ~ 0(1). Then, the fast character of the variable 
X is clarified, see similar arguments in [ >], and the system reads as 

\ at 

Our model reduction assumption consists in assuming that the curve defined 
by equation f{x,y) = is a good approximation when e ^ 1 to the slow 
manifold. This manifold coincides with the unstable manifold that joins the 
spontaneous point v^'^ to the two other stable equilibrium points {Si and 
5*3) after the bifurcation point while is part of the stable manifold before the 
pitchfork bifurcation. 

Due to the non-linearity of the function /, see (13) and (4), we cannot 
expect an explicit formula for /(x, y) = 0. Nevertheless, since dxf{0, 0) 7^ 0, 
the resolution in the neighborhood of the origin is given by the implicit 
function theorem. Hence we can define a curve: 

x = x*{y) (14) 

such that f{x*{y),y) = in a neighbourhood of the origin. We also note that, 
by construction the approximated slow manifold x*{y), implicitly defined by 
(14), intersects the exact slow manifold at all equilibrium points, i.e. where 
both / and g vanish (nullclines). Finally, we can conclude the slow-fast 
ansatz, replacing the complete dynamics by the one on the approximated 
slow manifold, and obtain the reduced one dimensional differential equation: 



y = g{x*{y),y)- 



3.2 Stochastic term 



We consider now the stochastic terms of system (1). When changing the 
variable form i/ to X, also the standard deviation of the considered Brownian 



10 



motion should be modified. Indeed the new variables x and y are linear 
combination of vi and V2. For instance, consider two stochastic differential 
eqautions: dvi = f3id(,i, where are two independent normalized white noises 
and Pi are the two standard deviations, and take a linear combination of vi 
and f2 with real constant coefficients ai,a2- x = aivi + 022^2- Then x must 
obey to the following stochastic differential equation: 

dx = V(ai/3i)2 + {a2(32Ydi. 

In our case, X = P~^{i' — v^'^), then we have: 

y = {P-%i{ui - vl'^) + {P'%2{1^2 - 4') 

or developing and considering dvi = /3id(^i, 

dy = {P-%iMii + {P'%2Mi2 

Since in our model /3i = /32 = /?, and considering the above discussion, we 
can write for a white noise : 

dy = P^{{P~^)2i? + {{P-^)22fdi. 
Finally, we conclude that the reduced one dimensional model reads: 

y = g{x*{y),y) + pydi (15) 

with j3y = j3\J {{P~'^)2if' + {{P~^)22)^ ■ We note that in the unbiased case 
j3y = (3, since P is given by (9). 



3.3 One dimensional Fokker-Planck model 

We can now consider the Fokker-Planck (or progressive Kolmogorov) equa- 
tion associated to the one dimensional stochatic differential equation (15). 
This gives the reduced dynamics on the approximated slow manifold x = 
x*{y). Let us remark that this reduced SDE is obtained except at the bi- 
furcation point and therefore valid whenever the slow-fast decomposition is 
verified or in other words whenever e is small. 

Consider the probability distribution function q{t,y), for t >0 and y G 
[~ym,+ym], then it must obey to the following one dimensional Fokker- 
Planck equation: 

dtq + dy I g{x*{y), y)q - -^dyQ 1 = (16) 
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with no-flux boundary conditions on y = ±ym- 

0^ 

g{x*{y),y)q- -^dyq = 0. 

Since equation (16) is one dimensional, it is always possible to find the 
effective potential G{y) being the derivative of the flux term g(x* {y),y). In 
other words, we can always define the potential function: 




Figure 5: Comparison of the potential —G for various values of AA = 0{top — 
left),Om{top - right), Om{bottom - left),0.1{bottom - right). 

Moreover, we can explicitly obtain the stationary solutions of (16), i.e. 
the solutions qs{y) independent on time t, as follows: 

'Zs(2/) = ^exp(-2G(y)//32), (17) 

with Z a suitable normalization constant. As explained also in [■]], these 
stationary solutions are the asymptotic equilibrium states for the solution of 
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the Fokker-Planck equation. In other words, letting time to go to infinity, 
the solution q{t,y) to (16) must converge to qs{y)- We have shown in [-3] 
that the decay to equilibrium for the two dimensional problem was exponen- 
tial. Nevertheless, this decay is so slow due to the small positive eigenvalue 
associated to the spontaneous state that the simulation shows metastable 
behavior for large times. Hence it is relevant to have a simple approximated 
computation of their asymptotic behavior without need to solve the whole 
2D Fokker-Planck equation which is provided by this effective ID potential. 

4 ID model vs. 2D model 




Figure 6: Comparison of the marginals in the new variable y, for dif- 
ferents values of the biasing paramter : AA = 0{top — left),0.01{top — 
right), 0.05{bottom — left),0.1(bottom — right). Blue line: the marginal 
computed by means of the ID problem. Black line: marginal computed 
from the 2D problem. Red line: the stationary marginal for the ID. Final 
time is 400 seconds. 

In this section, we numerically compare the solutions obtained for the one 
dimensional reduced Fokker-Planck equation (16) to the one of the original 
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two dimensional model (5). Concerning the numerical scheme for the two 
dimensional problem, we refer the reader to the detailed description in [ ]. In 
particular, we are interested in the solutions at equilibrium. As announced 
in section 3.3, we have an explicit formula for the solution at equilibrium in 
ID (17) by computing the primitive G{y). On the contrary, in the 2D set- 
ting, we cannot have such formulae and the computational time to approach 
equilibrium is very large. 





Figure 7: Comparison of the marginals in the variable vi for different values 
of the biasing parameter: AA = 0{top— left), 0.01{top— right), 0.05{bottom — 
left),0.1{bottom — right). Blue line: the marginal computed for the ID 
problem at equilibrium. Black line: marginal computed from the 2D problem 
for time T = 400s. 



In Figure 6, we plot the solution at equilibrium of the ID problem (the 
blue line) and compare it with the projection of the two dimensional solution 
on the new variable y (the black line). We remark that the black line is not 
too smooth since we are projecting a 2D distribution on a uniform quad- 
rangular mesh onto an inclined straight line. We can see a good matching 
in the unbiased case (AA = 0). In the biased cases, the results are differ- 
ent: for AA = 0.01, one clearly sees that even if we have computed until 
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the final time of 400 seconds, both the 2D and the ID solutions have not 
reached equilibrium and the 2D results are closer to equilibrium; while for 
AA = 0.05, or 0.1, the difference is smaller since the drift is strong enough 
to push all particles toward only one of the equilibrium points and there is 
only one population bump at least for the 2D results. The 2D results are 
closer to equilibrium at AA = 0.05 while at AA = 0.1 the ID are closer. 

On the other way round, we can also compare the marginals obtained 
from the two dimensional problem with the projections of the solution for 
the one dimensional problem on the vi and/or 1^2 axes. In figure 7 we show 
the comparison for various AA. Note that AA = 0.01 is the most interesting 
case as discussed in the previous figure. In fact, for larger AA, at equilibrium, 
the particles are almost all concentrated around one of the two stable points. 
Thus, no bump is visible around the second one (even in the one dimensional 
reduced solution), and for the unbiased case the matching is excellent. We 
warn the reader in order to compare Figures 6 and 7 that increasing values 
of y correspond to decreasing values of vi. 

The results demonstrate that the ID reduction is worth to obtain in- 
formation about the behavior at equilibrium. In the next section, we shall 
investigate the time dependent solution q{y,t) of equation (16). 

4.1 Time dependent solution 




Figure 8: Evolution in time of the distribution in the y variable, for the 
biasing parameter AA = 0.01 with snapshot every dt= 20 or 0.2 s (left) and 
20^ or 200s (right). 

We are here interested also on the time behavior of the solution q{y, t) of the 
ID Fokker-Planck equation (16). For instance, we may compare the time 
evolution of momentum for the 2D and the reduced ID problem. Thus, we 
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need to compute not only the stationary solution of equation (16), but also 
its time dependent solution. We choose to discretize equation (16) using 
implicit in time finite difference method. The evolution of the ID reduced 
model illustrates again the slow-fast character of this problem. In fact, we 
observe in Figure 8 the evolution in time of the density q{y, t) for small (left) 
and for large (right) times respectively. The convergence toward the final 
stationary state is quite slow compared to the fast division toward the two 
bump distribution at the initial stages. 

We describe now how to recover all the moments of the partial distri- 
bution function in the (2^1,^2) plane, using the probability distribution 
function q{y) solution of (16) and the approximated slow manifold x*{y). 

The function p^ is concentrated along the the curve v = (z^i(y), J^2(y)) 
given by v'^'^ + P{x*{y),y)^ , see (11). We parametrized this curve by means 
of a curvilinear coordinate and define 

V{y) = \\Pix'*iy),lf\\. 

Then, for any test function ^ = ^(z^i, z/2), the moment of the probability 
distribution function p^ = Pe{^i^^2) is defined by 

Mm = I ^pedvidu2 , 
Jo. 

and given by 

Mvt := j ■^{{vi{y),V2{v))q{y)dy. 

This formulae can be used to compute either classical moments of Pf, or 
marginals by choosing e.g. ^ = to get the f2-marginal as a function 

of //. Note that q has to be normalized in such way that its total mass (along 
the slow manifold) is equal to 1 i.e. Mij,=i = 1. 

Let us illustrate this metastability by the evolution of the first moments 
of the distribution in Figure 9. The initial data is a Dirac measure located 
above the spontaneous point (x = 0, y > small). We choose AA = 0, 
(5 = 0.3 and a;+ = 2.35. We use a implicit scheme in order to have no 
stability constraint on the time step. The number of discretization point 
is 200 and the time step is At = 0.01 for the left plot. It shows the fast 
dynamics: first the Dirac measure diffuses onto a Gaussian blob and moves 
quickly toward the spontaneous (unstable) state, then the Gaussian blob 
splits in the two bumps around the two stable equilibrium points. It seems 
that the solution has reached a equilibrium but it evolves very slowly. The 
figure on the right corresponds to At = 100 and shows this slow evolution 
toward the real equilibrium state. We will comment more below. 
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Figure 9: Evolution in time of the first order moments for the ID Fokker- 
Planck reduction. Left Figure: final time is 10^. Right Figure: final time is 
10^. Time unit is 0.01s 

We can finally compare the marginals in ui for the ID reduced model 
and the 2D simulations in Figure 10. We can conclude that the transients 
of the 2D are captured extremely well by the ID reduced model. 




Figure 10: Comparison of the marginals in variable vi, for the biasing pa- 
rameter A A = 0.01. Left Figure: ID reduced model. Right Figure: 2D 
simulation. 

5 Reaction time and Performance 

In the previous sections, we have numerically studied the accuracy of the 
reduced ID model with respect to the 2D original one. We discuss now some 
other information we can obtain from the ID problem, namely: the escaping 
time, section 5.1, and the probability density to belong to a sub-domain of 
the phase space, section 5.2. 
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5.1 Escaping time 



Fixed a bias AA and for a variable (3 we can easily compute the escaping 
time. In fact we recall that the Kramers formula [14]: 

E{t) ~ eMHc/f^^) 
where Hq is the maximum difference of the potential G 

Hg = Graax — Gmini (18) 

apply in the one dimensional framework, without needing to compute the 
solution q{t,y) of the Fokker-Planck equation (16). Recall that Gmax cor- 
responds to the potential value at the spontaneous state while Gmin corre- 
sponds to the minimum of the potential value at the two decision states. In 
Figure 11, we plot the potential gap computed by means of (18) as a function 
of the bias AA. 




Figure 11: Computed potential gap as a function of A A. 

In the 2D problem, since the drift is not the gradient of a potential, 
Kramer's rule does not apply and the escaping times can only be computed 
for the unbiased case, AA = 0. In fact, for AA = the problem is sym- 
metric in vi and 1^2 and thus, we know that the firing rates will separate in 
two identical bumps. Then, starting the computation from an initial data 
narrowly concentrated around one stable equilibrium point (say ^i), the es- 
caping time T can be defined as the time needed to have half of the total 
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mass moving to a neighborhood of S3. In particular, the expectation E(T) 
has an exponential behavior and its associated potential gap is Hg = 0.1. 
Of course, this kind of argument cannot be extended to the biased case and 
thus, the ID reduced model is essential. 

5.2 Probability densities - Performance 

We can compare the value of the probability pi of the firing rates to be 
in some domain for the 2D Fokker-Plank model and the ID reduced 
FP model. In particular, we shall compare: pp the probability for the 2D 
problem that at time t the firing rates belong to {{x,y) : y > 0}; with p^ 
representing the probability for the solution of the ID problem to belong to 
17+ = [0,ym]) see Figure 4. We fix the standard derivation /3 = 0.1, and let 
the bias AA varying from to 0.05, since the values of pp and p^ for AA 
bigger than 0.03 are already very close (the relative error being of the order 
of 10-4). 

We recall that, in the 2D problem, we has to wait for a very long time 
in order to reach equilibrium, since we have a meta-stable situation, see 
[3]. Nevertheless, we note that the Pp{t) profile is exponentially increasing 
converging to an asymptotic value poo- We can extrapolate this value from 
the values of Pp{t) for some initial iterations as follows. Assume that the 
probability pp behaves like: 

p{t) = Poo - aexp(-i/r) 

where p^o, a and r have to be determined by an "exponential regression". 
For a sequence of time ti (of the form tj = tQ + T * i, i = - ■ ■ N that 
corresponds to the computed values of pp), we define Api as the difference 
Pp{ti) - pp{ti + T), we get: 

Api = aexp(-fi/r)(l - exp(-r/T)). (19) 

Taking the log and the difference between two indexes i and j we obtain the 
expected value of r as: 

log(Ap,)-log(Apj) 
r ~ . 

Finally from (19), knowing r (and T), we recover a, and the asymptotic limit 
Poo is uniquely determined by: 

Poo = p{to) + aexp(-to/r). (20) 
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We show in Figure 12 the comparison between the values for the one di- 
mensional computation (red line) and the one extrapolated from the 2D 
computation, see equation (20), using a final time of 20 seconds (blue line). 
Note that the non-smoothness of the blue line (2D extrapolation) may be 
due to the fact that for computing pp we need to compute the inner product 
for any point v of the phase space (i^i, V2) ■ < (i^eq — i^), Pi,j >, and we choose 
for different values of AA the same equilibrium point I'eq and matrix P: for 
instance, for AA = 0.016, 0.018, 0.020, 0.022, 0.024 we choose the values of 
AA = 0.020: 

_ / 3.0448158 \ / 0.7003255 -0.6959201 \ 

^"'^ ~ \ 3.2397474 J ' ~ \^ 0.7138236 0.7181192 J ' 




0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 



Figure 12: Comparison of the values for p+ and poo with respect to the 
biasing parameter AA G [0, 0.05]. Red line: the values computed from the ID 
reduced problem. Blue line: the values extrapolated from the 2D problem. 
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